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Magnetic systems are an exciting realm of study that is being explored on smaller and smaller scales. One 
extremely interesting magnetic state that has gained momentum in recent years is the skyrmionic state. It is 
characterized by a vortex where the edge magnetic moments point opposite to the core. Although skyrmions 
have many possible realizations, in practice, creating them in a lab is a difficult task to accomplish. In this work, 
new methods for skyrmion generation and customization are suggested. Skyrmionic behavior was numerically 
observed in minimally customized simulations of spheres, hemisphere, ellipsoids, and hemi-ellipsoids, for typ- 
ical Cobalt parameters, in a range from approximately 40 nm to 120 nm in diameter simply by applying a 
field. 



I. INTRODUCTION 

A skyrmion, theorized first by Skyrme in 1962^, is a state 
with a vectorial order parameter which is aligned at the sys- 
tem boundary at an opposite direction to what the order pa- 
rameter assumes at the origin. Skyrmions may appear in di- 
verse arenas, such as elementary particles 1 5 , liquid crystals 6 , 
Bose-Einstein condensates^, thin magnetic filmPH quantum 
Hall systemjM Hl, a nd potentially vortex lattices in type II 
superconductor^^. Being able to experimentally observe 
or generate skyrmions is a current research thrusPGJD 

In this work we demonstrate via micromagnetic simulations 
that achieving a skyrmion is as simple as creating a nanopar- 
ticle of many possible geometries, which is large enough to 
support a single vortex but small enough to prevent multiple 
vortices. The demagnetization energy allows for the forma- 
tion of a vortex at zero-field. We find that as the field increases 
such that it lies in a direction opposite to the core, the mag- 
netization at the edges may realign itself parallel to the field 
direction more readily than the magnetization next to the core. 
Immediately prior to annihilation of the vortex (i.e., the flip- 
ping of the magnetization at the system core to become paral- 
lel to the applied field direction), the skyrmionic state is most 
notable. We observed this, relatively ubiquitous, effect in sys- 
tems with disparate geometries- spheres, hemispheres, ellip- 
soids, and hemi-ellipsoids. It may be possible to generalize 
this process so as to experimentally synthesize a skyrmion lat- 
tice by simply creating an array of nanoparticles with tunable 
size and spacing, such as by self-organzation 2021 . Prelimi- 
nary simulations of a two-by-two grid of Cobalt hemispheres 
of radius 20 nm with varying inter-hemisphere separation in- 
dicate that beyond a threshold distance of twice the radius, an 
array of skyrmions is formed. As the center to center separa- 
tion is steadily increased, the skyrmionic state becomes more 
lucid. For small separations, interactions partially thwart the 
creation of the individual skyrmions. 

As is well known, we can quantify a skyrmionic state by 



calculating the Pontryagin index (also known as a winding 
number) that is given by^ 
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In this expression, is the two dimensional anti- 
symmetric tensor and M is the normalized magnetization. 
For a single skyrmion, this winding number (or topological 
charge) is equal to unity. Skyrmions are characterized by the 
non-trivial homotopy class ^(S 2 ). This homotopy class is 
characterized by an integer that, for this case, is the Pontrya- 
gin index. States with different integer skyrmion number (the 
Pontryagin index) cannot be continuously deformed into one 
another. 

In the current context, the skyrmionic state resides on a two 
dimensional plane. On each spatial point of the plane, there is 
a three dimensional order parameter which, in our case, is the 
magnetization M. Topologically, a skyrmion is a magnetic 
state such that when it is mapped onto a sphere (via stereo- 
graphic projection) resembles a monopole or hairy ball. This 
means that on mapping from a flat space to the surface of a 
sphere, the individual magnetic moments will always point 
perpendicular to the surface of the sphere, much like a mag- 
netic monopole. 

The above topological classification is valid for an "ideal" 
skyrmion on an infinite two-dimensional plane or disk with 
the condition that the local moment M(r) at spatial infinity 
(irrespective of the location r on the infinite disk) all orient in 
the same direction: lim r ^ oc M(r) = M . In such a case M 
corresponds to the magnetization at the "point at infinity". On 
applying a stereographic projection of the infinite plane onto 
a unit sphere, Mo maps onto the magnetization at the north 
pole of the unit sphere while the oppositely oriented M at the 
origin corresponds to the magnetization at the south pole. In 
such a case, the winding number is identically equal (in abso- 
lute value) to unity. In many physically pertinent geometries, 
including the systems simulated in this work, there are finite 
size limits which only allow the magnetization M to exhibit 
the trend of approaching a uniform value Mq as one moves 
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away from the center of the system. In this case, the integral 
in Eq. [T] is not an integer. However, it is clear that, in the 
limit of infinite planar size, these states would become ideal 
skyrmions and the winding number Q would approach an in- 
teger value. 

The remainder of this article is organized as follows. In 
Section [TTJ we provide necessary background. We briefly de- 
scribe the simulations employed in this work and discuss ener- 
getic considerations. Section [III] reports on our central result- 
the numerical observation of skyrmions. We discuss a higher 
dimensional generalization and the possibility of generating 
skyrmion lattices. We conclude in section [TV] with a summary 
of our results. 



II. THEORY 

A. Simulation Theory 

In this work of simulating magnetic states of nanoparticles, 
the Object Oriented Micromagnetic Framework (OOMMF) 
1.2a distribution as provided from NIST was utilized^. The 
OOMMF code numerically solves the Landau-Lifshitz Ordi- 
nary Differential Equation given by, 



dM 
di 



-\^\MxH eff -^Mx^MxH eff ) (2) 

where M is the magnetization, 7 is the Landau-Lifshitz gy- 
romagnetic ratio, M s is the saturation magnetization, a is the 
damping coefficient, and H e ff is the effective field given by 
derivatives of the Gibbs free energy. The Gibbs free energy, 
in this case, is given by^, 



G = / ( ^ C [(^) 2+ (^) 2+ (^ 7 )^ 

~M-H' -M-H )dr (3) 

where a, f3, and 7 are the directional cosines, C is propor- 
tional to the exchange stiffness constant and depends on the 
crystal structure, w a is the crystalline anisotropy term, H' is 
the demagnetization field, and Ho is the external magnetic 
field. The crystalline anisotropy term can be expressed in 
terms of anisotropy constants, K\ and K2, and directional 
cosines as, 



w a = K x (a 2 (3 2 + /3 2 7 2 + 7 V) + K 2 a 2 (3 2 j 2 . (4) 

In the simulations, a metastable state was determined to 
have been reached when the maximum torque experienced by 
any one magnetic moment, measured in deg ^ es , dropped be- 
low 0.2. Once this level of torque was reached, the magnetic 
state data were saved to a file along with the other proper- 
ties of the system, including but not limited to, the energies 
associated with each contribution, overall magnetization, and 
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TABLE I. Table of parameters used in the simulations of particles 
in this work. The exchange stiffness constant, saturation magnetiza- 
tion, and crystalline anisotropy constant are material specific and are 
chosen for Cobalt. The damping constant, Landau-Lifshitz-Giblert 
gyromagnetic ratio, and stopping torque are material independent pa- 
rameters. 



number of iterations. The magnetic field was then changed 
to the next value and the iterations continued until saturation 
of the magnetization was obtained. The magnetic field steps 
were chosen such that half the steps (typically, a few hundred) 
were during the increasing field portion and the other half in 
the decreasing field portion. The data stored in the file were 
used later to generate the hysteresis plots, track the energy 
changes associated with the field variations, and the spatial 
orientations of the magnetic moments. Unless specified oth- 
erwise, the parameters chosen in the simulations correspond 
to those for Cobalt, as shown in Table |H 



B. Energy Considerations 

In our simulations, we considered field, demagnetization, 
and exchange energies. For simplicity, we neglected crys- 
talline anisotropy effects. The field tries to align the local 
magnetic moments parallel to it while exchange effects fa- 
vor an alignment of the magnetic moments with their near- 
est neighbors. The (universally geometry borne) demagneti- 
zation energy directly relates to dipole-dipole interaction^. 
Demagnetization energy is often the dominant term for long 
range behaviors while exchange effects tend to dominate at 
short spatial scales. 

As is well known, the competition between the long range 
and the short range energy contributions leads to the creation 
of domain walls. The demagnetization favors oppositely ori- 
ented moments at the expense of exchange effects that favor 
slow variations amongst neighbors. Ultimately, this tradeoff 
gives rise to domain walls in micromagnetic systems. 

The potential energy from demagnetization of a system is 
given by 
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where h\ is the effective field at position i that originates from 
all other dipoles. This field can be written as 



h' i =H' + -irM + tif, 



(6) 



where H' is the megascopic field from the poles due to M 
outside of a physically small sphere around site i. The second 
term subtracts the effective field inside an arbitrary small re- 
gion (or sphere) centered about point i, and h" is the field at 
site i created by dipoles inside this region. In general, h" de- 
pends on the crystal lattice structure. In the continuum limit, 
the sum becomes an integral of the form, 



M)dV. (7) 



The tensor A in the third term depends only on the crystal 
structure and local magnetization and can grouped with crys- 
talline anisotropy. This tensor also vanishes for cubic crystals 
identically. The second term in this expression is a constant 
proportional to M 2 and can be ignored. The A tensor also 
vanishes for cubic crystals identically leaving, 
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The demagnetization field, H', can equivalently be derived 
from Maxwell's equations. It can be expressed as the negative 
gradient of a potential, U that satisfies the equations, 



V 2 U in = lB V-M 

v 2 u out = 0, 

with the surface boundary conditions, 
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where the constant 7# is, in our units, 4ir. 

Lastly, the potential needs to be regular at infinity, such that 
\rU\ and \r 2 U\ are bounded as r — » oo. Our simulations 
directly capture the demagnetization field effects. 

From the standpoint of energy, for a skyrmion to be pos- 
sible, the dimensions of the ellipsoid must be larger than the 
critical dimensions at which vortices can nucleate in a given 
system. For example, for the hemispherical geometry, with 
the typical values of Table [TJ the critical radius was found to 
be 19 nm. For larger radii, vortices are the preferred state be- 
fore reaching zero field. The vortex will nucleate such that 
the core is parallel to the field and the remainder of the vor- 
tex lies in the plane perpendicular to the field. Once the field 
begins to oppose the direction of the moments at the core, the 
energy cost of eliminating the core is significantly higher than 



allowing the outer magnetic moments to align more with the 
field. When the exchange energy cost of the skyrmionic state 
becomes greater than the demagnetization energy for a uni- 
form magnetization, the core flips, annihilating the skyrmion, 
and the magnetization saturates. Immediately, prior to this, 
though, a skyrmionic state can be achieved. 

Ezawa 25 raised the specter of a skyrmionic state in thin 
films via the computation of the energy of such assumed varia- 
tional states within a field theoretic framework of a non-linear 
sigma model. Dipole-dipole interactions may stabilize such a 
state below a critical field. Our exact numerical calculations 
for the evolution of the magnetic states demonstrate that not 
only are skyrmionic states viable structures, but are actually 
the precise lowest energy state for slices of hemispheres and 
other general structures. 



III. RESULTS AND DISCUSSION 
A. Observation of a Skyrmion 

As our numerical simulations vividly illustrate, just prior to 
the annihilation of the vortex, the magnetic moments at the 
edge of the system start to orient themselves in a direction 
opposite to that in the core. On increasing the radius of the 
simulated hemispheres and spheres, the configurations next to 
the basal plane better conformed to the full skyrmion topology 
(i.e., that on an infinite plane).t should be noted here, that as 
the radius of a hemisphere increases, the crossover to a double 
vortex state will eventually occur, but if one vortex is main- 
tained, in the limit of large radii, a full skyrmion would be 
expected. This may be possible in materials with large ex- 
change constant and small saturation magnetization. In what 
follows, we will employ the typical values appearing in Table 
[I] The skyrmion state for the bottom layer (basal plane) of a 
hemisphere of radius 24 nm is shown in Figure [T] 




FIG. 1 . Vector plot of the skyrmion state for the bottom slice of a 
hemisphere of radius 24 nm. Not all local magnetic moments are 
shown for the sake of clarity. 

A similar configuration was observed in simulation runs for 
nanospheres. For a sphere, symmetry does not favor any par- 
ticular direction, but that symmetry is broken once a field is 
applied. Skyrmions were observed in runs of spheres large 
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enough to support a vortex which corresponds to a radius of 
« 15 nm. As the radius of the sphere increases, the edge mag- 
netic moments and the core magnetic moments become more 
antiparallel. A skyrmion in a sphere of radius 59nm is shown 
in Figure [2] 





FIG. 4. Vector plot of the skyrmion state in a hemi-ellipsoid with 
major axis of 20nra and minor axis of 15 nm. The slice is along the 
base of the hemi-ellipsoid. Only a subset of local magnetic moments 
is shown for clarity. 



0.8 T. These data are shown in Figure [5] 



FIG. 2. Vector plot of the skyrmion state in a sphere of radius 59nra. 
The slice is along the equator of the sphere. Only a subset of local 
magnetic moments is shown for clarity. 

Once skyrmions were observed in these systems, it begged 
the question, "Do these occur in ellipsoids and hemi- 
ellipsoids?" Upon examining this, indeed skyrmions can be 
observed in oblate ellipsoids and hemi-ellipsoids as shown in 
Figures [3] and [4] 




FIG. 3. Vector plot of the skyrmion state in an ellipsoid with major 
axis of 20nra and minor axis of 15nra. The slice is along the equator 
of the ellipsoid. Only a subset of local magnetic moments is shown 
for clarity. 

To verify that these are structures approach those of 
skyrmions and to quantitively monitor their deviations from 
an ideal skyrmionic state (for which the Pontryfin index is 
unity),we computed the Pontryagin index at different cross 
sections of the hemisphere. These cross sections were those 
of the hemisphere with planes parallel to the basal plane(i.e., 
that at the base of the hemisphere). For a hemisphere with ra- 
dius 30 nm, we calculated the skyrmion number Q for thirty 
individual parallel layers vertically separated by lnm. We nu- 
merically evaluated the integral of Eq. [T]for all of these layers 
and examined how it changes as the field increases from to 
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FIG. 5. Plot of the Pontryagin index versus the z-coordinate of the 
slice taken from the hemisphere of radius 30 nm. These are shown 
for increasing field from zero field (dark blue dot-dash line), 0.2 T 
(green dotted line), 0.4T (red dashed line), and 0.6T (teal solid line). 

Visualizing this in the geometry of the hemisphere specifi- 
cally, one can look at how the Pontryagin index varies along 
various planes of a hemisphere, starting from the equator and 
moving to the pole. It can be clearly seen that the skyrmionic 
behavior exists for most of the height of the hemisphere and 
only the cap deviates from the rest of the system. The size 
of this cap depends on the given field strength as can be seen 
in the case of field (Figure [6(a)] ) and with a field of 0.6 T 
(Figure 6(b) ). At higher fields, prior to the annihilation of the 
vortex, the Pontryagin index approaches an integer value, as 
expected for an ideal skyrmionic state. 

Performing similar analysis on the hemi-ellipsoids and vi- 
sualizing the Pontryagin index and its variance with height, it 
can be seen that the same behavior exists in a less extreme way 
than the hemispheres. This behavior can be seen in Figure [7] 
for hemi-ellipsoids of fixed 30 nm major axis and varying mi- 
nor axis. 



(a) 



(b) 



FIG. 6. Three dimensional plots of the Pontryagin index for a hemisphere of radius 30 nm at (a) zero field and (b) 0.6 T 



In examining the hysteresis behavior of the hemi-ellipsoids, 
one can see a trend as the z-dimension goes from the hemi- 
sphere radius (20 nm) to the minimum simulated size of 
5 nm. This trend shows a movement from extensive vortex 
and skyrmionic behavior in the more hemispheric geometries 
and less vortex and skyrmionic behavior in the more ellip- 
soidal geometries. 

Although it will not be considered in this work, crystalline 
anisotropy could influence the formation of skyrmions in a 
number of ways. In the case of a single crystal, the vortex state 
would be more difficult to nucleate and thus the skyrmionic 
state is less energetically favorable. When many crystalline 
grains are present, the results discussed here are valid as the 
large number or randomly oriented crystals will, on average, 
not favor any direction, and thus will not favor any one direc- 
tion. 



B. Generalization to a Hedgehog 

These results lead to the question of whether this can be 
generalized to more than two dimensions. The natural gen- 
eralization from the two-dimensional skyrmion to a three- 
dimensional magnetic state would be the hedgehog. The 
hedgehog resides in three spatial dimensions coupled with a 
three dimensional order parameter. The canonical example 
of a hedgehog isM = M s r where the magnetization always 
points outwards. A skyrmion is related to a hedgehog via a 
stereographic projection from the sphere onto a plane where 
the south pole of the hedgehog projects to the core of the 
skyrmion on the plane and the north pole of the hedgehog 
projects to the points at infinity on the plane. Calculating the 
demagnetizing field for this state in a sphere gives rise to a 
potential and field equal to 

U{r)= lB M s {r-R), (13) 
H = -j B M s r. (14) 



Plugging this into Equation [81 one finds the energy of the 
hedgehog to be 2nM^ (An/ 3) Br. Comparing this to the en- 
ergy of the uniformly magnetized state, (1/2) (47r/3) 2 M^R 3 , 
it can easily be seen that the hedgehog has three times the en- 
ergy of the uniform state. This, combined with the fact that 
the exchange energy and the field energy will favor the uni- 
form state, the hedgehog state will not be possible in a sphere. 

If one were to continuously deform the hedgehog by ro- 
tating the local magnetic moments by tt/2 such that M = 
M s f(z)(j) where f(z) is a function that goes to as z —> 
such that the exchange energy does not diverge, one would 
find the demagnetization energy of that state to be identically 
0. The field energy in this system is also for a field that 
is applied along the z-axis. The exchange energy is given by 
(4tt/3)RC where C is the exchange stiffness constant. The 
total energy of this state is equal to the exchange energy, and 
comparing this to the uniform state, a hedgehog of this form 
is favorable for, 

For C = 2.5 x 10" 11 J/m and M s = 1.4 x 10 6 A/m as it 
is for Cobalt, at field, this radius works out to be « 3.5/im. 



C. Skyrmion Array 

It is illuminating to consider the possibility of an array of 
skyrmions. As briefly discussed below, we find that effective 
particle interactions may thwart the creation of a skyrmion lat- 
tice when these particles are not far separated. However, for 
sufficiently large center to center separations, a Skryme lat- 
tice may be achieved. In preliminary simulations of arrays of 
nanoparticle arrays, simulations of a two-by-two grid of hemi- 
spheres of radius 20 nm with a variable separation show that 
a center to center separation of four times the radius is close 
enough that the nanoparticles still interact magnetically and 
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FIG. 7. Plots of the Pontryagin index and how it varies with height inside hemi-ellipsoids of 30 nm radius major axis as the minor axis varies 
from 15 nm (a) to 10 nm (c) to 5 nm (e). This is shown for (a) field equal to 0.2 T pointing in the negative z-direction (perpendicular to the 
face of the hemi-ellipsoids). As will be noted, the existence of skyrmionic behavior is not prevalent in the more flattened hemiellipsoinds and 
vanishes at this field between minor axis 15 nm and 10 nm. The associated partial hysteresis loops for each of these hemi-ellipsoid runs are 
shown in Figs, (b), (d), and (f), respectively. 



prevent the formation of an array of skyrmions. As expected, 
further separation should approach the the single particle re- 
sult of skyrmions, as we briefly discuss next. 

The transition from the array of particles which support in- 
dividual vortices to the array of particles that are clearly in- 
teracting with each other can be seen in Figure [5] In this fig- 



ure, the annihilation of the vortices can be seen as the par- 
ticles realign their magnetization to form a state where the 
local magnetization orients in the counterclockwise direction 
from particle to particle, yet within each particle, when mov- 
ing in the counterclockwise direction, the local magnetization 
changes from oriented in the negative z-direction to the posi- 
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tive z-direction. 

In repeating these simulations for a 3x3 array of hemispher- 
ical nanoparticles, the same behavior was observed. This ar- 
ray was similar to the 2x2 array in that it had nanoparticles 
with diameters of 40 nm and center to center separation of 
80 nm. The annihilation of the vortices occurred at a slightly 
larger field (0.08T rather than 0.1T) as shown in Fig|9] 



IV. CONCLUSION 

We conclude with a brief synopsis of our findings. We 
carried a systematic numerical study of the magnetization of 
small nanoparticles in the presence of an external magnetic 
field. These systems were simulated for different sizes and 
geometry (sphere, hemisphere, ellipsoids). Our analysis ig- 
nored anisotropy (crystalline, shape, strain, etc.) effects. We 
find that, as has been widely reported in the literature^^J be- 



yond a critical diameter, the particles enter into a single vortex 
state under zero external field; multiple vortices are possible 
for much larger particles. Our key new result concerns the 
creation of skyrmions in the single vortex state. As the field 
is increased, vortex annihilation is accompanied by the for- 
mation of a skyrmionic state wherein the magnetization of the 
vortex core points to a direction opposite to that at the edge 
of the nanoparticle. Our result illustrates how geometry plays 
a pivotal role. Spheres and hemispheres more readily achieve 
skyrmionic states than higher eccentricity ellipsoids. Our pre- 
liminary results suggested that for center to center separations 
larger than twice the particle diameters, an array of skyrmions 
may be realized. More detailed studies of skyrmion lattices 
for such particle arrays are planned for the future. 
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FIG. 8. Vector plot of a 2x2 array of hemispheres with radius 20 nm and center to center separation 80 nm at fields of 0.12 T pointing in 
the negative z-direction (a) and 0.1 T pointing in the negative z-direction (b). Colorscale corresponds to z-component of the local magnetic 
moment in units of A/m. 




FIG. 9. Vector plot of a 3x3 array of hemispheres with radius 20 nm and center to center separation 80 nm at fields of 0.1 T pointing in 
the negative z-direction (a) and 0.08 T pointing in the negative z-direction (b). Colorscale corresponds to z-component of the local magnetic 
moment in units of A/m. 



